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ABSTRACT 

The standard procedure to generate initial conditions in numerical simulations of structure 
formation is to use the Zel'dovich approximation (ZA). Although the ZA correctly reproduces 
the linear growing modes of density and velocity perturbations, non-linear growth is inaccurately 
represented, particularly for velocity perturbations because of the ZA failure to conserve momen- 
tum. This implies that it takes time for the actual dynamics to establish the correct statistical 
properties of density and velocity fields. 

We extend the standard formulation of non-linear perturbation theory (PT) to include tran- 
sients as non-linear excitations of decaying modes caused by the initial conditions. These new 
non-linear solutions interpolate between the initial conditions and the late-time solutions given 
by the exact non-linear dynamics. To quantify the magnitude of transients, we focus on higher- 
order statistics of the density contrast S and velocity divergence O, characterized by the Sp and 
Tp parameters. These describe the non-Gaussianity of the probability distribution through its 
connected moments {5P)c = Sp( ( 6^ )c = Tp ( 6^ f-\ We calculate Sp{a) and Tp{a) to 

leading-order in PT with top-hat smoothing as a function of scale factor a. 

We find that the time-scale of transients is determined, at a given order p, by the effective 
spectral index rtoff. The skewness factor ^3 (T3) attains 10% accuracy only after a « 6 (a w 15) 
for UeS ~ 0, whereas higher (lower) rics demands more (less) expansion away from the initial 
conditions. These requirements become much more stringent as p increases, always showing slower 
decay of transients for Tp than Sp. For models with density parameter f2 7^ 1, the conditions 
above apply to the linear growth factor, thus an fl = 0.3 open model requires roughly a factor 
of two larger expansion than a critical density model to reduce transients by the same amount. 
The predicted transients in Sp are in good agreement with numerical simulations. 

More accurate initial conditions can be achieved by using second-order Lagrangian PT (2LPT), 
which reproduces growing modes up to second-order and thus eliminates transients in the skew- 
ness parameters. We show that for p > 3 this scheme can reduce the required expansion by more 
than an order of magnitude compared to the ZA. Setting up 2LPT initial conditions only requires 
minimal, inexpensive changes to ZA codes. We suggest simple steps for its implementation. 



Subject headings: cosmology: theory; cosmology: large-scale structure of universe; methods: 
analytical; methods: numerical 
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1. Introduction 

Gravitational instability is widely considered responsible for the formation and evolution of large-scale 
structures in the Universe. The non-linear nature of gravitational clustering makes analytical calculations 
only possible for models with restrictive symmetries or in cases where the density fluctuations are small 
enough that a perturbative approach is possible. For this reason, numerical simulations are a vital resource 
for understanding how large-scale structures form and evolve in the Universe. 

The standard procedure in numerical simulations is to set up the initial perturbations, assumed to be 
Gaussian, by using the Zel'dovich (1970) approximation (Klypin & Shandarin 1983; see also Efstathiou et 
al. 1985, hereafter EDFW). This gives a useful prescription to perturb the positions of particles from some 
regular pattern (commonly a grid or a "glass" ) and assign them velocities according to the growing mode in 
linear perturbation theory. In this way, one can generate fluctuations with any desired power spectrum and 
then numerically evolve them forward in time to the present epoch. 

Although the Zel'dovich approximation (hereafter ZA) correctly reproduces the linear growing modes of 
density and velocity perturbations, non-linear correlations are known to be inaccurate when compared to the 
exact dynamics (Grinstein & Wise 1987, Juszkiewicz, Bouchet & Colombi 1993, Bernardeau 1994, Catelan 
& Moscardini 1994, Juszkiewicz et al. 1995). This implies that it may take a non-negligible amount of 
time for the exact dynamics to establish the correct statistical properties of density and velocity fields. This 
transient behavior affects in greater extent statistical quantities which are sensitive to phase correlations 
of density and velocity fields; by contrast, the two-point function, variance, and power spectrum of density 
fluctuations at large scales can be described by linear perturbation theory, and are thus unaffected by the 
incorrect higher-order correlations imposed by the initial conditions. 

In this work, we therefore concentrate on higher-order statistics such as the one-point cumulants of 
density and velocity divergence fields, characterized by the so-called Sp and Tp parameters (Goroff et al. 1986; 
see Eq. ( p7|) below for their definition). We use non- linear perturbation theory (PT) in order to provide a 
quantitative description of transients. For this purpose, we extend the standard formulation of PT to include 
the full time-dependence of density and velocity solutions to arbitrary non-linear order in the perturbation 
expansion. We show that initial conditions can be thought as exciting non-linear decaying modes in the 
evolution of perturbations which lead to transient behavior. Although we assume the initial conditions given 
by the ZA, the present formalism can be useful to explore other non-Gaussian initial conditions as well. We 
present results up to p = 8, although the techniques developed in this work make calculations possible to 
higher orders if desired with minimal extra effort. Given that current angular surveys are able to measure 
up to 5*9 (Gaztahaga 1994; Szapudi, Meiksin & Nichol 1996), and that the situation in redshift surveys will 
greatly improve in the near future, it is important to address the issue of what requirements are needed in 
order to determine accurately these statistical measures of clustering in numerical simulations. 

The problem of transients from ZA initial conditions is well known in the PT literature, and it has 
been pointed out many times before (e.g., Juszkiewicz, Bouchet & Colombi 1993, Juszkiewicz et al. 1995). 
However, the only quantitative understanding of the magnitude of this problem so far comes from the 
numerical work of Baugh, Gaztafiaga & Efstathiou (1995, BGE hereafter), who studied the transients in 
the skewness of the density field and found that for CDM simulations, an expansion a « 3 away from the 
initial conditions was necessary to erase the memory of the ZA-induced skewness. Recently, the effect of 
transients in ^3 at large scales was seen in the numerical renormalization solution by Couchman & Peebles 
(1997, Figure 12), where fluctuations at large scales are applied by the ZA at each renormalization step and 
then expanded by a factor of two in scale factor. For higher order moments and velocity fields moments, 
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there is no analysis (numerical or otherwise) available in the literature. In this sense, this work constitutes a 
natural development of the subject and should be useful in designing numerical simulations of cosmological 
structure formation. 

We now give a somewhat detailed summary of the contents of this paper for those readers not familiar 
with the technical machinery of non-linear PT and mostly interested in the relevance of the results for 
numerical simulations. 



In Section 2.1 we review the standard formulation of non-linear PT. Eqs. (|9D-(p_ 1\) show how to obtain 
non-linear solutions to the equations of motion by a recursive procedure from linear PT solutions, as first 
derived by Goroff et al. (1986). 



Section 2.2 presents the description of ZA initial conditions in terms of Lagrangian trajectories (as used 



to set initial conditions in numerical simulations, Eqs. (p^-(p^)) and in Eulerian space (Eqs. (14)-(|6| 
which is in fact the most convenient way to calculate the statistical properties of initial conditions. We 
explore two different ways of setting initial velocities, the standard ZA approach which leads to Eq. (|l] 
and the EDFW scheme which yields Eq. (fsh for the divergence of the initial peculiar velocities. 



In Section 2.3, we extend the non-linear PT formalism to include a description of transients. The main 
result in this regard is the recursion relation in Eq. (|2^) which gives the non-linear solutions to the equations 
of motion including the properties of initial conditions (represented by the first term) and how they excite 
non-linear decaying modes (transients). These solutions interpolate between the initial conditions (ZA) and 
the late-time exact-dynamics recursion relations given by Eq. (^t|). 

In Section |^, explicit expressions for the Sp and Tp are given in terms of the solutions presented in the 
previous section. A simple way to understand the physical meaning of the Sp parameters is provided by 
the Edgeworth expansion (e.g. Juszkiewicz et al. 1995; Bernardeau & Kofman 1995) for the probability 
distribution function (PDF) of the standardized variable i' = S/a 

P{v) = PoM {l + ls, H,{v) a H,{v) + ^ Hd^)] ^2 + . . . }, (1) 

where PQiv) = (27r)^^/^ exp(— 1/^/2) is the Gaussian distribution, and Hn{v) are the Hermite polynomials, 
e.g. H^v) = —iv. An equivalent expression can be written down for the velocity divergence PDF involving 
the Tp parameters. We see from Eq. (|l|) that as cr —> 0, the PDF approaches a Gaussian distribution. The 
first-order correction to Gaussianity is given by the skewness factor which, being proportional to the third 
moment of the PDF (see Eq. (p7|)), measures the asymmetry between overdense and underdense regions. 
Gravitational clustering leads to a positive Sz (Peebles 1980), which means that the high-density tail of the 
PDF is enhanced and the underdense tail is suppressed with respect to a Gaussian distribution, as expected 
from the attractive nature of the gravitational force. Similarly, since underdense regions expand and occupy 
a larger fraction of the volume, the most probable value of 6 becomes negative, from Eq. (|l|) the maximum 
of the PDF is in fact reached at 

<^max ~ --^ o-, (2) 

to first order in a. We thus see that the skewness factor contains very useful information on the shape of 
the PDF. Similarly, higher-order Sp parameters provide further information on the PDF shape, for example, 
the next term in Eq. (|^) contains a contribution from the kurtosis factor 5*4 that gives the lowest-order 
measure of flattening of the PDF tails relative to a Gaussian distribution. Comparison with exact PDF 
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results when available (e.g. for the ZA where the PDF can be calculated non-perturbatively, see Bernardeau 
& Kofman 1995) and numerical simulations (Juszkiewicz et al. 1995) shows that a few terms in Eq. (|^) are 
good enough to reproduce the shape of the PDF near its maximum. On the other hand, the Edgeworth 
expansion rapidly breaks down at the tails of the distribution, even becoming non-positive definite, because 
they receive appreciable contributions from the whole hierarchy of Sp, and the Edgeworth series cannot be 
truncated at finite order. In fact, the use of powerful generating function techniques (Balian & Schaeffer 
1989; Bernardeau 1992, 1994; Bernardeau & Kofman 1995), which can be thought as a resummation of 
the Edgeworth expansion, shows that the tails of the PDF in the asymptotic region are determined by the 
infinite series of Sp parameters through their generating function. 

In summary, the main point of this digression is to convince the reader of the physical relevance of the 
Sp and Tp parameters as statistical measures of clustering. It should be clear from this discussion that an 
accurate determination of these is desired in numerical simulations in order to represent the statistics of 
large-scale structures correctly (BGE; Gaztafiaga & Baugh 1995; Lokas et al. 1995; Juszkiewicz et al. 1995; 
Colombi, Bouchet & Hernquist 1996; Bernardeau et al. 1997). Regarding the calculation of Sp and Tp from 
leading-order (tree-level) PT (that is, at lowest non- vanishing order in the variance of the density field), 
it turns out that vertices (defined in Eq. (|3^)) contain all the necessary information about the non-linear 
dynamics. A recursion relation for vertices that includes transients is presented in Appendix 

In Section^ we present the perturbative results for transients in Sp and Tp. Figure^ shows the magnitude 
of transients for 5^ (3 < p < 8) for different spectra, and similarly for Tp in Fig. ^. A comparison with 
numerical simulations with initial velocities as in EDFW is presented in Fig. |[ Figures |l^ and ^present the 
corresponding transients in Sp and Tp when initial conditions are generated using second-order Lagrangian 
PT (2LPT). Note the significant improvement when compared to the equivalent plots in Figs. |l| and || for 



Section |^ contains our conclusions and a final discussion. Technical material regarding details of the 
calculations and further analytic results are presented in Appendices |^, and Finally, Appendix |^ 
presents a simple procedure to implement second-order Lagrangian PT initial conditions in numerical simu- 
lations which requires minimal numerically inexpensive modifications to a standard ZA code. 



In the single-stream approximation, prior to orbit crossing, one can adopt a fluid description of the 
cosniological iV-body problem, where the relevant equations of motion correspond to conservation of mass 
and momentum and the Poisson equation (e.g., Peebles 1980). Since vorticity decays in an expanding 
universe, the system can be described completely in terms of the density contrast S(x) = [p(x,t) — p\/p 
and the velocity divergence, 9 = V ■ v, where v(x) is the peculiar velocity. Defining the conformal time 
T = J dt/a, where a(r) is the cosmic scale factor, and the conformal expansion rate Ti. = dhia/dr = Ha, 
where H is the Hubble constant, the equations of motion in Fourier space become (e.g. Scoccimarro & 
Frieman 1996) 



the ZA. 



2. Dynamics 



2.1. Standard Formulation of Perturbation Theory 




(3a) 
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^^^^ + h{t) 0{k,T) + ^nn^{T) ~6{k,T)^~ j d'ki j d'k^ [&]2 l3{k,kiM) o{ki,T) 9{k2,T), 

(3b) 

where we used the Fourier transform representation S{k) = (27r)^'^ J cPxS{x) e~^^'^. In Eqs. (^), k denotes 
a comoving wave number, [i5d]2 = 6-D{k — ki — k2), with the Dirac delta distribution, and 

a(fe, fei) ^ /3(fc, fci, fc2) ^ (4) 

which describe mode-couphng due to the non-hnear dynamics. Equations (^) and (^) are vahd in an 
arbitrary homogeneous and isotropic universe, which evolves according to the standard Friedmann equations. 



In linear PT, the solution to the equations of motion ( |3a| ) and (3b) is given by 

S{k,T) - D{T)6,{k), (5) 
0{k,T) = -H{t) f{n,A) D{t) 5i{k), (6) 

where D{t) is linear growing mode, which from the equations of motion must satisfy 
and f{n,A) is defined as 

a In a n dr 

ExpHcit expressions or D{t) and /(fi. A) are not needed for our purposes (see e.g. Peebles 1980), but useful 
fits in the cosmologically interesting cases are f{n, A) « rj3/5 when A = (Peebles 1976), and f{n, A) « f^^/g 
when n + r^A = 1 (Bouchet et al. 1995). When f2 = 1, the linear growth factor D{t) becomes the scale factor 
a(T) and / = 1. In this case, the PT solutions at each order n scale as a"(r) and a general recursion relation 
is available that gives the PT solutions at arbitrary order (Goroff et al. 1986, see Eqs. ( |ll| ) below). When 
0.^1 and/or A 0, the PT solutions at each order become increasingly more complicated, due to the fact 
that growing modes at order n in PT do not scale as a"(T). Furthermore, the solutions at each order become 
non-separable functions of t and fc (Bouchet et al. 1992, 1995; Bernardeau 1994b; Catelan et al. 1995), 
and there appear to be no general recursion relations for the PT kernels in an arbitrary FRW cosmology. 
However, a simple approximation to the equations of motion, fi//^ — 1, noted by Martel & Freudling (1991) 
for second order PT, leads to separable solutions to arbitrary order in PT in which the linear growth factor 
D{t) plays the role of the scale factor, and the same recursion relations as in the Einstein-de Sitter case 
are obtained for the PT kernels (Scoccimarro et al. 1998). All the information on the dependence of the 
PT solutions on the cosmological parameters ft and A is then encoded in the linear growth factor, D{t), 
which in turn corresponds to the normalization of the linear power spectrum. Equivalently, in this case the 
equations of motion can be written independently of cosmology by taking the linear growth factor as a time 
variable using Eq. (||) (Nusser & Colberg 1997). 

The PT solutions for arbitrary cosmological models in this approximation can be written down as 
(Scoccimarro et al. 1998, Appendix B.3) 



~6{k,T) = ^i^"(T) 5„(fe), (9a) 

oc 

d{k,T) = -H{t) f{n,A)J2D"{T) d^{k), (9b) 
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The equations of motion then determine 5n{k) and On{k) in terms of the Hncar fluctuations, 



5„(fc) = y d3fci...y d3fc„ [fc]„Fi^)(fci,...,fe„) '5i(fci)---,5i(fe„), (lOa) 
en{k) = f d^h... f d^kn[Su]nGi:\ki,...,kn) 5i{ki)...S,{kn). (10b) 



where [S-o]n = [S-d] (A; — fei — ■ • • — fe„), and the fimctions -F'i*'' and Gn ^ are constructed from the fundamental 
mode coupHng functions a{k,ki) and P{k, ki, k2) by a recursive procedure (Goroff et al. 1986; Jain & 
Bertschinger 1994), 



n-l 



Fn{ki,. ■ ■ ,k„) = ^ 



Gn{k 



1 , . . . , «,„ J 



Gm{kl J • ■ ■ ; ^Tjx) 



, (2n+l) «(fc,fc(™)) i^„_,„(fc„,+i,...,fc„) 

(2n + 3)(n - 1) 

+2 /5(fc, fc^ \ fe^ -') G„-m(A;„i_|_i, . . . , fc„) 

Gm{kl , • • • , fcm) 

(2n + 3)(n- 1) 

-\-2n [}(k^k^ \k^ ^) Gn — m {km+l j ■ ■ ■ i ^n) 



E 



3 Q;(fc, fe^ ^) Fn—mikm+lj ■ ■ ■ j 



(11a) 



(lib) 



(where fc^"^ = fci + . . . + A;„, fc^""") ee fc™+i + . . . + fc„. A; = fe^"' +A;("""\ and Fi = Gi = 1). Symmetrizing 
Fn and G„ over the wavectors ki, . . . ,k„ leads to Fn^ and gI^^ as needed in Eq. (p^). 



2.2. Initial Conditions: Zel'dovich Approximation 

Let us now discuss the perturbative formulation of the Zel'dovich (1970) approximation, which we shall 
use to calculate the statistical properties of the initial conditions. In the ZA, the motion of each particle 
x{q, t) at initial position q is obtained by applying linear PT to its Lagrangian displacement, '4'(<7, r), 

x{q,T)^q+^{q,T)^q-D{T) V^'^'Hq), (12) 

where 0^^^ (q) denotes a Lagrangian potential given by the initial conditions, see Appendix ^ for a more 
detailed discussion of ZA in Lagrangian space. This implies that the velocities of particles initially at q are 
given by 

-Z?(t) W(t) / V0(i)(g). (13) 

In Eulerian space, it can be shown that the ZA is equivalent to the dynamics obtained by replacing the 
Poisson equation by the ansatz (Munshi & Starobinsky 1994, Hui & Bertschinger 1996) 

v(x, t) = '^^—V^ix, r), (14) 

which is the relation between velocity v{x,t) and gravitational potential ^{x,t) valid in linear PT theory. 
The ZA therefore fails to conserve momentum, i.e. the Euler equation will only be satisfied to linear order. 
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When calculating one-point cumulants of density and velocity fields, is convenient to use the perturbative 
version of ZA, since the non-perturbative probability distribution function (PDF) is singular due to caustic 
formation in the sense that one-point cumulants ( 5^ )c of the density field diverge for p > 1 due to the high- 
density tail (e.g., see Bernardeau & Kofman 1995). The perturbative formulation of ZA follows from replacing 
Eq. (|l^) into the equations of motion (^, which yields the corresponding recursion relations (Scoccimarro 
& Fricman 1996) 



i^„2A(fci,...,fc„) 



E 



m— 1 



n(n — 1) 



(n-1) a(fc,fc(")) ^^„2f^(fe™+i,...,fc„) 
+P{k, k^"^\ '"■') G^^„(fem+i, . . . , kn) 



(15a) 



^ n— 1 

G^^(fci,...,fc„) = VGf„A(fci,...,fc„)/3(fc,fc(™),fc("-'"))G^^„(fc,„+i,...,fe„). (15b) 

Upon symmetrization these give a simple result for the density field PT kernels (Grinstein & Wise 1987) 



, hfi) — 



1 (fe • fei) (fc • fc„) 



h.2 



(16) 



In numerical simulations, the ZA is generally used to set the initial perturbations as follows (see e.g. EDFW). 
A Gaussian density field, 5i{k) is generated in Fourier space from the desired linear power spectrum, and 
therefore the Lagrangian potential (f)'^^^{k) ~ —Se{k)/k'^ is Fourier transformed and its gradient taken to 
yield V q(l)^^\q) at each "unperturbed" particle position, denoted by q, which is usually described by a grid 
or a "glass". This gives the displacement field which is used in Equation (|l^) to displace the particles 
from their unperturbed positions and imposes the ZA density field from the initial Gaussian field, i.e. 
— M[5i{q)], where M denotes the ZA mapping imphed by Eqs. (10a) and (p^). Equation (^ can 



then be used to assign the velocities to particles. In this case, the velocities then satisfy 

d{x) = -fn f^D'^T) j d'ki... j <fkn [SD]nGl^{ki,...,kn) (5,(fci) . . . (5,(fe„). 

n=l 

It is in fact straightforward but tedious to show that Eq. ( p^ in Lagrangian space implies Eq. (|l^) in Eulerian 
space with the kernels G^^ given by Eq. (|l5b[). 



(17) 



Although most existing initial conditions codes use this prescription to set up their ZA initial conditions, 
there is another prescription to set initial velocities suggested by EDFW, which avoids the high initial 
velocities that result from Eq. ( p^ ) because of small-scale density fluctuations approaching unity when 
starting a simulation at late epochs (low redshift). This procedure corresponds to recalculate the velocities 
from the gravitational potential due to the perturbed particle positions, obtained by solving again Poisson 
equation after particles have been displaced according to Eq. (^2|) . Linear PT is then applied to the density 
field to obtain the velocities, which implies instead of Eq. (^ and Eq. (|l^) that 



9{x) 



-fU Y.D-{r) 

n=l 



(Fki 



d'K [fc]„ Ff^^{k^, . . . , fe„) 5,(fci) . . . 5t{kn). (18) 
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Therefore, in this case, the initial velocity field is such that the divergence field Q(x) = 9(x)/{—f TL) has the 
same higher- order correlations as the ZA density perturbations. Comparing Eq. ( p^ and (|l^) we see that 
both prescriptions agree in Unear PT, where F^^ — = 1, but they differ at second and higher-order PT, 
since / G^^ for n > 1. This imphes that these two different alternatives of setting the initial velocities 
affect the magnitude of the transients, as we shall see below. It turns out that in fact, the prescription given 
byEq. 



is closer to the exact dynamics given by Eqs. (10b) and (lib) than Eq. (|l7|), therefore, the second 
method will excite less decaying modes and consequently transients will decay faster than in the standard 
ZA scheme. 



2.3. Transients in Perturbation Theory 

Perturbation theory describes the non-linear dynamics as a collection of linear waves, (5i(fc), interacting 
through the mode-coupling functions a and (i in Eq. (^. Even if the initial conditions are set in the growing 
mode, after scattering due to non-linear interactions waves do not remain purely in the growing mode. In the 



standard treatment, presented in Section p.l| , the subdominant time-dependencies that necessarily appear 
due to this process have been neglected, i.e., only the fastest growing mode (proportional to D") is taken 
into account at each order n in PT. In this subsection we generalize the previous results to include the full 
time dependence of the solutions at every order in PT. This is necessary to properly address the problem of 
transients from ZA initial conditions. 

In this case we write the perturbative solutions as (see Eqs. (|^)): 



l{k,r) = ;^n"(T) 5„(fe,i?), (19a) 

oo 

0>, r) = -U{t) f{n. A) e(fe, r) = -H{t) f{n, A) J2 D'\t) 0^{k, D), (19b) 

71- 

where Sn and 9n are written in terms of PT kernels as in Equation 1^) 



n=l 



$i")(fc,D) = / d^h... / d^kn [Mn ^l")(fci,...,fc„;i?) Si{ki)---5,{K), (20) 



where a = 1,2 and <I>^" = 6n, $2 = ^n- "^^^ kernels Ta now depend on the linear growth factor D{t) 
and reduce to the ones in the previous section when transients die out, that is J-'f^'^ — > F„, J^j"'' ~^ 



when D{t) — ^ 00. Also, Eq. ( poD incorporates in a convenient way initial conditions, if we assume that at 
D = l, A'"^ — xj^\ we effectively reduce the non-Gaussianity of initial conditions to a Gaussian problem, i.e. 
where the linear solutions Si{k) are Gaussian random fields. The kernels xi"' describe the initial correlations 
imposed at the start of the simulation, for the standard ZA scheme we have (see Eq. (pT|)) 

j{n) ^ pZA^ j{n) ^ qZA^ (21) 



whereas for velocities set from perturbed particle positions we have instead (see Eq. (18)) 

j{n) ^ pZA^ j{n) ^ pZA_ (22) 
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The recursion relations for !F^\ which solve the non- linear dynamics at arbitrary order in PT, can be 
obtained by replacing Eq. ([2C|) into the equations of motion, which yields 



n-l „2 

^ •^i™'' ('''ii • ■ ■ J ^m! s) JF^ \km+iT ■ ■ ■ s), (23) 

where z = In Z?(r) and we have assumed the summation convention over repeated indices, which run between 
1 and 2. In Equation (03) the matrix 



gab{z) = — 



3 2 
3 2 



(24) 



is the linear propagator (Scoccimarro 1997b). The first term represents the propagation of linear grow- 
ing mode solutions, where the second corresponds to the decaying modes propagation. That is, in linear 
perturbation theory, density and velocity perturbations propagate in time according to 



$W(z)=5a.(z) 1>i'^(0). 



(25) 



If initial conditions are set in the growing mode, then <I>i^''(0) oc (1,1) vanishes upon contraction with the 



second term in gab{z), whereas the first term reduces to the familiar ^^a\z) = e^$a^''(0) = D{t) <I'q^''(0) 
The "scattering matrix" ^abc encodes the non-linear interactions and is given by 



7i2i(fe, fei,fc2) = a(fc,fci), (26a) 
l222{kMM) = l3{kMM). (26b) 



and zero otherwise, with k = ki + k2- Note that since gab{z) — > 5ab as z 0+, in Eq. ( p3| ) the kernels T^"^ 
reduce to Xa"'' at Z? = 1, where the initial conditions are set. Equation ( |2^ ) reduces to the standard recursion 
relations in Eq. ( p|) for Gaussian initial conditions (xl"^ = for n > 1) when transients are neglected, i.e. 
the time dependence of JFa"'' is neglected and the lower limit of integration is replaced by s = — oo. Also, 
it is easy to check from Eq. (2^) that if Xa"'' ~ {FmGn), then JT^"-* = {Fn,Gn), as it should be. In what 
follows, the kernels in Eq. (E3) are assumed to be symmetrized over its arguments. Note that these kernels 
are no longer a separable function of wave- vectors and time. 

Equation (^3|) gives useful insight into the nature of nonlinear PT solutions. For example, second order 
solutions are built from the interaction (represented by the 7 matrix) of two linear waves that are propagated 
freely to the present time using Eq. (^4|), plus the contribution from initial conditions propagated to the 
present represented by the first term in Eq. (p^). After scatterings, waves do not remain purely in their 
growing modes, and from Eq. (^3|) one can check that there is a contribution from decaying modes propagation 
even for the fastest growing mode included in standard treatments of PT. More importantly for the present 
purposes is the fact that the contribution to the PT solutions at a given order n and "time" z depend on 
all the n-scattering processes that happened between s = z and s — 1, where initial conditions are set. By 
assuming that initial conditions are set in the "infinite past" (s = —00), the standard formulation of PT 
presented in the previous section contains no information on the time scale that takes to erase the memory 
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of the correlations imposed by the initial conditions. A more detailed treatment of this formalism and other 
applications is left for a future paper (Scoccimarro 1997b). 



3. Statistics 

We are interested in one-point cumulant statistics for the density and velocity divergence fields, which 
can be characterized by the Sp and Tp parameters defined as 

where R is the smoothing scale, and the subscript c denotes the connected contribution, i.e. the gaussian 
value is subtracted off. Note that we define cumulants for the velocity divergence field in terms of 8 (a;) as 
defined in Eq. ( |l9b[ ), which differs from the standard convention by a factor (—1)^. In this work, we use 
top-hat smoothing, which is described by a window function W'Tii{kR) 

WMx) = ^ = A [sin(a;) - xcos(a;)] (28) 

in the Fourier domain, where Ju{x) is a Bessel function. To simplify the notation, we henceforth denote: 

daj EE 47r fcf dk, P{ki) W,,„j = WTnilh + ... + kj\R), (29) 

where P{k) is the power spectrum, and denotes the variance of the density field fiuctuations 

( S{k)5{k') )=SB{k + k') P{k) a^{R) = J d^k P{k) W^^^ikR). (30) 

A systematic framework for calculating correlations of cosmological fields in PT has been formulated 
using diagrammatic techniques (Fry 1984; Goroff et al. 1986; Wise 1988; Scoccimarro & Frieman 1996). From 
this point of view, leading order PT for the statistical quantities of interest corresponds to tree graphs, next- 
to-leading order PT contributions can be described in terms of one-loop graphs, etc. These diagrammatic 
techniques assume Gaussian initial conditions, although in principle they can be extended to any non- 
Gaussian model by adding the appropriate new vertices (Wise 1988). In this work we are interested in a 
very particular non-Gaussian initial condition, that given by the ZA as usually imposed to start numerical 
simulations. As shown in the previous section, in this case one can include the non-Gaussianity directly into 
the non-linear solutions and therefore the problem reduces effectively to one dealing with Gaussian initial 
conditions. 

We now write down the standard expressions for the Sp parameters in terms of non-linear kernels. In 
the discussion that follows, analogous analysis always applies to the velocity divergence. In the following we 
shall use the scale factor a(r) to denote the time dependence, hut in models with density parameter 17 7^ 1 the 
same equations will he valid upon replacing a{T) hy the linear growth factor D{t). For the skewness factor, 
we have 

S3{R,a) = J Pid^'ki P2d^k2 W1W12W2 T[^\ki,k2;a). (31) 
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The calculation of tree-level diagrams such as this and the expressions that follow is simplified by using the 
fact that tree-level PT corresponds to taking the spherical average of PT kernels (Bernardeau 1992, 1994b). 
We therefore define the angular averaged smoothed kernels wi"^ by (a = 1, 2): 



1 



— (i/„,/i„) = t^i' 



(«) 



/dill 



Wl...n ^fHfci,...,fc„), 



(32) 



n\ ' ' " J \ m J \ 47r 
whose recursion relations for top-hat smoothing are obtained directly from Eqs. ( p3|) in Appendix A. The 
vertices m„ and ^„ are the (Eulerian) smoothed counterpart of those defined in Bernardeau (1992). Recently, 
Lagrangian smoothed vertices have been defined by Fosalba & Gaztahaga (1997). The skewness factor 
smoothed at scale R can then be rewritten as: 



53 (i?, a) 



da I da 2 



3 vi{xi;a) V2{xi,X2;a) vi{x2;a), 



(33) 



where Xi = kiR and i^i(a;) = W{kR). To simplify the notation, we henceforth suppress the dependence on 
the scale factor a. Similarly, the final expression for the kurtosis factor is 



Si{R) 



da\ dal dal 



12 Vi{xi) y2{xi,X2) V2{X2,X3) Vi(x3) + A iyi{xi) Vi{x2) Vlixz) U3{xi,X2,X3) 

(34) 



whereas for Sc, (the "pentosis" parameter according to the nomenclature of Chodorowski & Bouchet 1996), 
we obtain 



S,{R) = 



dal da2 da^ da"^ '' 
%R) 



60 J^l(xi) 1^2{X1,X2) V2{X2,X3) V2{x3,X4) Vi{xi) + 60 Vi{xi) Vi{x2) 
X V'i{xi,X2,X'i) V2{X3,X4) Vi(xa) + ^ Vl{xi) Vi{x2) Vlixz) Vi{xa) Vi{xi, X2, x^, xa) , (35) 

and, finally, for 5*6 we have 



Sq[R) 



dal c?cr| dal ^'^4 ^^5 



6 J^i(xi) 1^1 (2:2) i^iixs) iyi{x4) iyi{x5) ly^ixi, X2, X3, xa, X5) 



a^°{R) 

+ 120 i^i(xi) 1^1(2:2) viixs) i'a{xi,X2,xz,xa) 1^2(2:4, 2:5) vi{xz) 

+ 90 Vi{xi) Vi{x2) 1^3(2:1, 0:2,2:3) V3{xz,XA,X,i) vi{xa) viix^) 
+ 360 Vi(xi) Vi{x2) l'z{xi,X2,Xz) V2{x3,Xa) 1'2{xa,X^) vi{xz) 

+ 360 vi{xi) V'i{xi,X2,X'i) V2{x2,xa) i^2(a;3,a;5) i^i(2;4) i^i{x5) 

+ 360 l^l{xi) V2{XI,X2) V2{X2,X^) V2{x3-,Xa) V2{xa,X<=,) Vi{xz) 



(36) 



The tree-like structure of these expressions is quite obvious. The combinatorial factors in each term 
corresponds to the number of labellings of each particular tree diagram; see e.g Fry (1984) for 3 < p < 6 
and Boschan, Szapudi & Szalay (1994) for up to p = 8. The expressions for Sr{R) and Ss{R) will not 
be reproduced here, but note that they can easily be obtained from the tree diagrams together with their 
combinatorial coefficients shown in table 1 of Boschan, Szapudi & Szalay (1994). The equivalent expressions 
for Tp are generated by simply replacing v by /i in the formulas above. 
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4. Results 

Despite the complicated appearance of the expressions given in the previous section, they can be cal- 
culated in a straightforward manner thanks to the special geometrical properties of top-hat smoothing 
(Bernardeau 1994, see also Hivon et al. 1995). Gaussian smoothing does not share these properties and will 
not be considered here (see Lokas et al. 1995). The recursion relations for the smoothed vertices Vn and /x„ 
as a function of scale factor a and smoothing scale R are derived in Appendix A from the recursion relations 
given by Eq. (p^). These vertices depend on scale R through derivatives of the window functions, which 
are then converted into derivatives of the variance with respect to scale, the 7p(i?) parameters defined by 
(Bernardeau 1994) 

lAR) - -^nA^- (37) 

Using the results in Appendix A, we get for the skewness parameters (7 = 71 = n^s + 3) 



„ / ^ Ir. 1 r34 1 1 / 26\ 12 34 6 12 , 

= a[^-^] + {y-^} + a(^-T) + 3^"-^-5^+35;?7I' « 

lr„ 1 r26 1 1/ 16\ 18 26 6 18 , 

where we have assumed ZA initial velocities, Eq. (|lj). On the other hand, for initial velocities set from 
perturbed particle positions, as in Eq. (|l8|), we have: 



^ / ^ Ir. 1 r34 1 1 / 22\ 16 34 2 16 , 

- a[^-^] + {y-^} + a(^-T)-3^"-^-5^-35^' ^''^^ 

rr^ f X Ip. 1 r26 1 1 / 22\ 24 26 2 24 

- a[^-^] + {y-^} + a(^-T) + 3^-y-^-5^+35^- ^^^^^ 

For r2 1, these expressions and the ones that follow are valid upon replacing the scale factor a by the linear 
growth factor D{t). For scale-free initial conditions, with spectral index n, the results for top-hat smoothing 
are restricted to n < 1, since for n = 1 the variance in top-hat spheres diverges and the logarithmic derivative 
7 becomes meaningless (the same restriction applies to p > 3 results). The first term in square brackets in 



Eqs. (38) and ( |39D represents the initial skewness given by the ZA (Bernardeau 1994), which decays with the 
expansion as (Fry & Scherrer 1994). The second and remaining terms in Eqs. ( ^ ) and ( p9| ) represent 
the asymptotic exact values (in between braces; Juszkiewicz, Bouchet & Colombi 1993, Bernardeau 1994) 
and the transient induced by the exact dynamics respectively; their sum vanishes at a = 1 where the only 
correlations are those imposed by the initial conditions. Similar results to these are obtained for higher- 
order moments, we refer the reader to Appendix B for explicit expressions. Note that for scale-free initial 
conditions, the transient contributions to Sp and Tp break self-similarity. 

Figure |l] shows these results as a function of scale factor a for different spectral indices, assuming that 
velocities are set as in the ZA. The plots show the ratio of Sp{a) to its "true" asymptotic value predicted 
by PT, 5'p(oo), for 3 < p < 8. The values at a = 1 correspond to the ratio of ZA to exact dynamics S'p's, 
which becomes smaller as either p or n increases. For the skewness, it takes as much as a = 6 for n = 
to achieve 10% of the asymptotic exact PT value, whereas spectra with more large-scale power, where the 
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ZA works better, require less expansion factors to yield the same accuracy. A similar constraint should hold 
for the bispectrum, the three-point function in Fourier space (Peebles 1980; Fry, Melott & Shandarin 1993; 
Scoccimarro et al. 1998). As p increases, however, the transients become worse and at p = 8 an expansion 
by a factor a = 40 is required ior n — to achieve 10% accuracy in Ss- This suggests that the tails of the 
PDF could be quite affected by transients from initial conditions. Furthermore, for models where < 1, 
the requirements on a translate into requirements on the linear growth factor D{t), which implies a more 
stringent constraint on a, i.e. an f2 < 1 simulation should be started earlier (at a higher redshift) than an 
f2 = 1 model. For example, an open model with Q = 0.3 typically requires a factor of two higher initial 
redshift than for = 1 (see Figs. ^ and || below). 

Figure ^ shows the corresponding results for the velocity divergence Tp parameters. We see that the 
effects of transients in this case is more severe than in the density field case, in particular as n increases, since 
the initial a = 1 values imposed by the initial conditions become quite different from the asymptotic exact 
PT values. For example, n — requires a « 15 for 10% accuracy in T3 (more than a factor of two larger 
than for S3). The situation quickly deteriorates as p increases. Again, the shape of the PDF of the velocity 
divergence should be quite sensitive to the presence of transients. Moreover, since statistics of the density 
field in redshift space contain contributions from velocities correlations, one expects that the redshift-space 
density field PDF will be more affected by transients than in real space. 

Figures ^ and ^ show the equivalent results for initial velocities set from perturbed particle positions, 
Eq. (p^). Comparing to the corresponding results in Figs. |l| and ||, these results show considerable improve- 
ment on the amount of expansion required to erase transients, by factors of two or three depending on the 
spectrum. We see that at most a « 3 is necessary to achieve 10% accuracy in 6*3, at least a factor two better 
than in the ZA velocities case. These results are in agreement with the numerical study of BGE for CDM 
models, with velocities assigned as in Eq. (^), in which they found that a w 3 was needed to recover the 
PT prediction for S3 at scales where rtcff ~ — 1. 

Figure |^ presents a comparison of the perturbative predictions for transients in Sp parameters with 
the standard CDM numerical simulations measurements of BGE, kindly provided by E. Gaztafiaga. The 
simulation evolved lOO'^ particles in a box 300 h~^ Mpc a side. Initial velocities are set according to Eq. ( [l^ ) 
as described in EDFW. The error bars in the measurements correspond to the variance over 10 realizations. 
We plot these N-body results by taking the ratio to the tree-level exact dynamics value predicted by PT, 
which has the advantage of reducing the main scale dependence. If there were no transients and no other 
sources of systematic uncertainties, all the curves would approach unity at large scales, where tree-level PT 
applies. Unfortunately, there are other sources of systematic uncertainties as we shall discuss below. 

The different symbols correspond to different outputs of the simulation: open triangles denote initial 
conditions (a = 1, ag = 0.24), sohd triangles (a = 1.66, erg — 0.40), open squares (a = 2.75, as = 0.66), 
and solid squares (a — 4.2, (Tg = 1-0). The results for S3 in the top left panel are those presented in 
Fig. 9 of BGE. For the initial conditions measurements (open triangles) there is some disagreement with 
the ZA predictions, especially at small scales, due to discreteness effects, which have not been corrected 
for. The initial particle arrangement is a grid, therefore the Poisson model commonly used to correct for 
discreteness is not necessarily a good approximation (sec BGE for further discussion of this point). The 
second output time (solid triangles) is perhaps the best for testing the predictions of transients: discreteness 
corrections become much smaller due to evolution away from the initial conditions, and the system has not 
yet evolved long enough so that finite volume corrections are important. For S3 we see excellent agreement 



with the predictions of Eq. (39a), with a small excess at small scales due to non-linear evolution away from 



the tree-level prediction. For p > 3 the numerical results show a similar behavior with increased deviation 
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at small scales due to non-linear evolution, as expected (see discussion below). For the last two outputs 
we see a further increase of non-linear effects at small scales, then a reasonable agreement (at least for ^3 
and 5*4) with the transients predictions, and lastly a decrease of the numerical results compared to the PT 
predictions at large scales due to finite volume effects, which increase with erg, R, and p (Colombi, Bouchet 
& Schaeffer 1994; BGE; Colombi, Bouchet & Hernquist 1996, Munshi et al. 1997). For reasons of clarity, 
only measurements with reasonable error bars have been included, essentially as in BGE (see their Fig. 13). 

The A^-body results show a systematic overestimate of the PT predictions at small scales, i.e. i? < 10 
h~^Mpc for the last two outputs. Here we should stress that the PT predictions in this paper correspond to 
tree-level quantities, i.e. valid in the limit of vanishing variance, cr^{R) —> 0. For a small variance, one can 
use one-loop PT to calculate the correction to the tree-level results, and neglecting transients one finds that 
in general (for n < — 1 in the scale-free case) 

SpiR) ^ Sl°^ +a^iR) Sl^^ + (40) 

where Sj^^ denotes the tree-level value and Sp'^^ the one-loop correction coefficient (Scoccimarro & Frieman 

1996) . For p = 3 and top-hat smoothing, = 3.86, S'f ^ = 3.18 for n = -2 (Scoccimarro 1997). For 

p = 4 there are no available PT results to one-loop, but in the spherical collapse approximation (that has 

been recently extended to take into account loop-corrections by Fosalba & Gaztahaga (1997)), 5'!°^ = 27.56, 
(2) 

^4 = 63.56 for neS = —2. This approximation neglects tidal effects, but comparison with numerical 

simulations and exact PT results shows a good agreement for the Sp parameters (Fosalba & Gaztanaga 

(2) 

1997) . For example, the skewness one-loop coefficient 5*3 = 3.21 for n = —2 is in excellent agreement with 
the exact PT result quoted above. These results show that at scales of a few h~^Mpc for CDM models, 
where rics ~ —2, one-loop corrections play a significant role even for cr^(i?) < 1, and increasingly with p 
(Fosalba & Gaztafiaga 1997). This is the reason for the excess of scale dependence at i? < 10 h~^Mpc in 
Fig. ^ Another interesting point of this exercise with one-loop corrections is to show how a late start in a 
simulation can mimic a "false agreement" with tree-level PT: transients tend to decrease the measured 5'p's, 
whereas one-loop corrections due to the finite variance tend to cancel this decrease. This may lead to the 
illusion that tree-level PT has a wider range of applicability than is really the case. 

In Figure]^, the PT calculations include the full dependence on 7p(i?) parameters. An approximation 
is sometimes used to calculate Sp for p > 3 in which ^p{R) — for p > 2, on the grounds that this is 
true for scale-free spectra and CDM models have a slowly varying spectral index. Figure ^ addresses the 
validity of this approximation for standard CDM models. The top left panel shows n^s = 7 — 3 and the 
^p{R) parameters for 2 < p < 5, whereas the top right panel shows the ratio of Sp parameters (4 < p < 7) 
calculated using ^p{R) = for p > 2 to the full calculation. We see that the approximation works quite well 
at small scales, but as R increases (and the regime where tree-level PT holds is reached) the approximation 
jpiR) — gradually breaks down. The bottom panels in Fig. ^ show the same calculation for a CDM model 
with shape parameter F ~ 0.21 (which just corresponds to a shift in scale with respect to the F — 0.5 model), 
and the corresponding calculation in the standard CDM model for the Tp parameters. In this latter case, 
the approximation is worse than for the 5*^ parameters. Note that we have suppressed the plotting of Tp 
ratios beyond the point where they become negative for reasons of clarity. 

The next two figures. Figs. ^ and ||, show the predictions of transients as a function of scale factor 
a for Sp and Tp parameters respectively, at smoothing scales R = 10, 100 h^^Mpc (left and right panels, 
respectively) for F = 0.5 and 0.21 CDM models (top and bottom panels, respectively). These calculations 
include the full dependence on jp{R) parameters, and correspond to ZA initial velocities. A similar set of 
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plots is presented in Figs. |9| and ^for the case of velocities set as in EDFW. Comparing these two set of plots 
we arrive to a similar conclusion than for scale-free spectra, velocities set from perturbed particle positions 
lead to a faster rate of convergence to the asymptotic exact values than just pure ZA initial velocities. The 
latter requires a factor of about two to three more expansion away from the initial conditions to reduce 
transients in density and velocity statistics by the same amount. We also see that for R = 100 h~^Mpc, the 
transients are more important than for i? = 10 h~^Mpc, as expected from the results in the scale- free case. 
In the bottom panels, we include upper horizontal axes which denote the scale factor a for a cosmological 
model in which il{a) = 0.3. This serves to illustrate the fact that for models with < 1, transients persist 
longer because the growth of fluctuations is governed by the linear growth factor which evolves more slowly 
than the scale factor. Therefore, to achieve the same accuracy regarding transient behavior, < 1 models 
should be started at a higher redshift than Q = 1 simulations. The results in Figs. 0-^show that an Q = 0.3 
model should be started at about a factor of two larger in redshift than anfl — 1 simulation. We remind the 
reader that this result is approximate; it depends on the assumption that f'^ /fl ~ 1, but this approximation 
is better than 25% for O > 0.3. 

Finally, in view of these results, one is led to ask whether it is possible to decrease the magnitude 
of transients besides the obvious solution of starting the simulation early enough. A natural candidate to 
improve upon the ZA to set the initial conditions is to use second-order Lagrangian FT (hereafter 2LPT, e.g. 
Moutarde et al. 1991, Buchert 1992, Bouchet et al. 1995). This procedure requires minimal additional com- 
putational cost over the standard ZA scheme, as discussed in detail in Appendix D. Since 2LPT reproduces 
growing modes of density and velocity perturbations to second-order, there are no transients in the evolution 
of S3 and T3 parameters in this approximation. Figures |ll| and ^ present the predictions for transients from 
2LPT initial conditions for different spectral indices as a function of scale factor a for 3 < p < 8. The details 
of this calculation and the analytic results are summarized in Appendix C. Compared to the corresponding 
results in Figs. |l| and ^ respectively (note the difference in scales), we see that 2LPT initial conditions lead 
to an improvement of more than an order of magnitude in the amount of expansion necessary to erase tran- 
sients over the standard ZA scheme. Moreover, it yields about a factor of four improvement with respect 
to ZA initial conditions with velocities set from perturbed particle positions; that is, a simulation started 
with 2LPT initial conditions at redshift Zstait = 10 would roughly correspond to a simulation started at 
-Zstart = 40 using the ZA with velocities as in EDFW, and ^start = 100 for a standard ZA procedure. Given 
that generating EDFW velocities and 2LPT initial conditions require similar additional computational costs 
over the Z A scheme (see Appendix ^) , in any case very small compared to the actual cost of running the 
simulation, 2LPT seems the best alternative to standard ZA methods. 



5. Conclusions and Discussion 

In this paper we give a perturbative analysis of the problem of transients from initial conditions when 
measuring moments of density and velocity fields in numerical simulations, where initial conditions are usually 
set using the Zel'dovich approximation (ZA). Although the ZA correctly reproduces the linear growing modes 
of density and velocity perturbations, non-linear growth is inaccurately represented by the ZA, because of 
the ZA failure to conserve momentum. This implies that it takes a non-negligible amount of time for the 
correct dynamics to establish the expected tree-level correlation hierarchy predicted by perturbation theory 
at large scales. 

We focussed on one-point cumulants of the density and velocity divergence fields, characterized by the 
Sp and Tp parameters. We extended the standard formulation of perturbation theory to include to arbitrary 



- 16 - 



order the transient behavior of non-hncar solutions that encode the information on the amount of time needed 
to overcome the correlations imprinted by the initial conditions. Using these results, we calculated the full 
time-evolution of Sp{a) and Tp{a) to tree-level with top-hat smoothing as a function of scale factor a. These 
results interpolate between the initial values set by the ZA and the asymptotic values expected from the exact 
dynamics at large scales. More importantly, they provide a quantitative estimation of the magnitude of the 
effect and the amount of expansion needed to achieve a given accuracy in the determination of moments of 
the density and velocity fields. Needless to say, there are many other uncertainties when measuring statistics 
such as one-point cumulants in numerical simulations which should be properly taken into account; this has 
been extensively discussed in the literature (e.g. Colombi, Bouchet & Schaeffer 1994, 1995; BGE; Colombi, 
Bouchet & Hernquist 1996; Szapudi & Colombi 1996; Munshi et al. 1997). 

We found that the magnitude of transients is determined, at a given order p, by the effective spectral 
index noff. If initial conditions are set at Oq = 1, obtaining the skewness factor S3 (T3) within 10% accuracy 
requires a ~ 6 (a ~ 15) for Tiofr ~ 0, whereas higher (lower) n^g demands more (less) expansion away 
from initial conditions. Furthermore, these requirements become much more stringent as p increases, always 
showing slower decay of transients for Tp than Sp, due to velocity correlations being poorly represented 
by the ZA. For models with density parameter < 1, the transients take more expansion factors to die 
out, since the relevant dynamical quantity that controls the growth of structure is the linear growth factor, 
which evolves more slowly than the scale factor. This implies, for example, that an open model where the 
final state corresponds to 17 = 0.3 requires roughly a factor of two larger expansion away from the initial 
conditions to erase transients than an 51 = 1 model. Thus, in general, numerical simulations of models with 
< 1 should be started at higher redshift than critical density models to reduce transients by the same 
amount. We have also explored the influence of setting initial velocities on the magnitude of transients, and 
found that velocities set as in EDFW from the gravitational potential due to perturbed particle positions 
(i.e., after the ZA displacement has been applied), the time-scale of transients is reduced by a factor of two 
or three depending on the spectrum. 

The results of the predicted transients in Sp for 3 < p < 6 were compared with measurements in standard 
CDM numerical simulations by Baugh, Gaztafiaga & Efstathiou (1995). We found good agreement, especially 
for intermediate output times where discreteness and finite volume effects are not important and tree-level 
PT applies over a wider range of scales. Our results show that the dependence of transients on spectral index 
is opposite to that of finite volume effects, which decrease with riog (e.g., Colombi, Bouchet & Hernquist 
1996; Munshi et al. 1997). In the comparison shown in Fig. ^ we see that at late times finite volume effects 
start to dominate over transients, systematically decreasing S5{R) and Sq{R) for i? > 20 h~^Mpc. However, 
as the size of a CDM simulation box is increased and higher spectral indices are probed, transients eventually 
dominate over finite volume effects, since the latter become quite small as — > 1. 

Another interesting issue regarding transients in numerical simulations is to investigate to what extent 
transients are present at smaller scales. As R decreases, the PT approach used here breaks down (as shown in 
Fig. ||for i? < 10 h~^Mpc in the last two outputs) since we use leading order (tree-level) PT, which yields the 
Sp and Tp parameters in the limit of vanishing variance. However, at least in the intermediate regime where 
the variance is of order unity and one-loop PT works well (Scoccimarro & Frieman 1996; Scoccimarro 1997), 
the corrections to the predictions in this paper take the form given by Eq. (|40|). In this case, the expected 
result is that transients will take longer to decay than for tree-level quantities because loop-corrections depend 
on higher-order non-linearities which are increasingly underestimated by the ZA. That is in fact the reason 
why transients in Sp and Tp take longer to decay as p increases. On the other hand, for CDM models the 
spectral index decreases at smaller scales, and the magnitude of transients will decrease. As smaller scales 
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are approached and the variance becomes larger than unity, shell-crossing becomes important and it is less 
clear what to expect, this depends on how much the dynamics of virializing high-density regions couples to 
the large-scale modes. 

Since the results of transients in this work concern the statistical properties of density and velocity fields, 
it is likely that other statistical measures of clustering will show a similar effect, in particular those sensitive 
to phase correlations, where non-linear dynamics provides the leading contribution. An obvious candidate 
closely related to the discussion in this paper is the measurement of one-point cumulants in redshift space 
(Hivon et al. 1995), which can be thought as appropriate cross-correlations between density and velocity 
fields (Scoccimarro, Couchman & Frieman 1997). Based on the present results, one expects in this case 
that the Sp parameters in redshift space will be more affected by transients that the corresponding ones 
in real space, since velocity correlations are more affected by transients. Other statistics sensitive to large- 
scale phase coherence, such as percolation studies, may show similar effects regarding transients to the ones 
investigated in this work. 

The effect of transients also implies that one must be very careful when comparing non-linear approx- 
imations and perturbative results to numerical simulations which have not been evolved for a long enough 
time. In the former case, comparison among non-linear approximations usually concludes that the ZA (or 
modifications thereof) is a good approximation to the fully non-linear numerical simulation. These conclu- 
sions should be taken with caution if the simulation has not been started early enough so that transients from 
the initial conditions are still present. Also, as discussed in Section ^ when comparing numerical simulations 
to PT predictions, the effect of transients combined with finite values of the variance characteristic of a late 
simulation start tends to create the false impression that tree-level PT has a wider range of validity than is 
really the case. 

The results in this paper should perhaps be viewed as a useful guideline for designing numerical sim- 
ulations with interest in measuring higher-order statistics of density and velocity fields, as well as other 
measures of clustering. The main conclusion in this regard is the requirement for an early enough start 
to avoid undesired effects from transients, particularly for studies of clustering statistics at high redshift. 
Although an early start does not cause significant overhead in the time required to run a simulation because 
early time-steps run quite rapidly due to weak clustering, there are reasons to avoid starting a simulation at 
very high redshift. The faster Hubble expansion rate at earlier times demands successively shorter time-steps, 
but more importantly, the numerical integration of the equations of motion leads to suppressed growth of 
small-scale modes that becomes more significant by starting at higher redshifts (Couchman 1997). For these 
reasons, it is desirable to use a better approximation to generate initial conditions. A natural candidate is 
second-order Lagrangian PT (2LPT), which reproduces growing modes to second-order in non-linear PT. 
As shown in this work, this can reduce the expansion required to decay away transients by more than one 
order of magnitude compared to the standard ZA scheme. This is a significant improvement, particularly 
given that its implementation is simple and only requires minimal inexpensive modifications to widely used 
ZA initial conditions codes. 

This work was in part motivated by a conversation with Hugh Couchman and P.J.E. Peebles. I would 
also like to thank Stephane Colombi, Pablo Fosalba, Josh Frieman, Dmitry Pogosyan, and Istvan Szapudi 
for comments and discussions, and Enrique Gaztahaga for providing numerical simulation measurements and 
for comments as well. Additional thanks are due to Hugh Couchman for numerous helpful discussions . I 
thank the Aspen Center for Physics for hospitality during the workshop "Precision Measures of Large-Scale 
Structure" , where this project was started. 
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A. Tree-Level Recursion Relations for Smoothed Vertices including Transients 

We define the angular averaged smoothed kernels wi"'' : 

and similarly for the angular averaged initial conditions kernels ryi"'' upon replacing J-l by . As shown 
in detail by Bernardeau (1994), the operation of smoothing one-point cumulants with top-hat windows can be 
easily calculated due to special geometric properties. In terms of the scattering matrix elements in Eq. ([2^), 
the following properties of angular integration hold 



(A2a) 



-W1W2 + -W1 k2R Wi + -W2 kiR W[, (A2b) 
3 6 6 



where the prime denotes a derivative (Bernardeau 1994). To prove these results, all we need is the expansion 
of top- hat windows such as W12, in terms of Wi, W2 (and their derivatives), and Legendre polynomials of 
the angle between fci and fe2- This, it turns out, is a straightforward application of Gegenbauer's addition 
theorem for Bessel functions (e.g., Watson 1944). Then, using Eq. ( A2) and the recursion relations in Eq. ( ^ ) 
we find: 



(n) / \ —nz I \ (n) \ •' f , n(s — z) I ^ / \ To (™) (n—m) (m) o (n—m) 

[z) = e gab[z) % ' + / ds e ^ U - gai{z - s) 3 ' + ^^2 di^^i 



+ g ga2{z - S) 



O , ,('") , ,("-'") I , ,(™) a , ,("-™) I Pi , ,(™) , ,("-™) 
Z UJ2 ^2 ' ^2 ^K^2 ^K^2 ^2 



(A3) 



where k = Ini? and the quantities in square brackets are evaluated at time s. This recursion relation (and 
its derivatives) can be used to obtain the smoothed vertices to arbitrary order in PT in terms of top-hat 
window functions and their derivatives with respect to scale, which are then straightforwardly converted to 



derivatives of the variance with respect to scale, yielding the 7p(i?) parameters defined in Eq. (37). If we 



ignore transients, and assume Gaussian initial conditions, Eq. (A3) reduces to (n > 1) 



, ,(") — J /X f'T,') {"i , /™) / ,("-™) I , ,(™) , (n~m) 



m— 1 



+ ^ (Ja2{n) 




(A4) 



with 



o-abin) 



(2n + 3)(n- 1) 



2n+ 1 
3 



2 

2n 



(A5) 
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If we further assume no smoothing, Eq. (A4) yields 



n-l 



(A6) 



where 7121 = 1, 7222 = 1/3, and zero otherwise, which is the hmit of Eq. (A2) as the smoothing scale 
goes to zero, R ~> 0. This simple recursion relation is an alternative method to reproduce the well-known 
results for unsmoothed vertices usually obtained using the spherical collapse approximation, i.e. 1^2 = 34/21, 
^3 = 682/189, 1^4 = 446440/43659, . . . , ^2 = 26/21, ^3 = 142/63, fi4 = 236872/43659, and so on (Bernardeau 
1992). 



B. Results for Transients in Kurtosis and Pentosis Factors 

In this appendix, we present results for the transients in 6*4, T4, 6*5 and Tc, from ZA initial conditions, for 
standard ZA velocities and initial velocities as in EDFW, set from perturbed particle positions (Eq. (|l8|)). 
For p = 4, evaluation of Eq. (HJ) using the methods described in Appendix A leads to 



60712 62 7 7 72 2 72 816 28 7 184 1312 8 7 1504 192 
64(a) = ' , , , 



1323 3 3 3 35 a 5a 75 a2 245 ai 5ai 4725 at 1225 a^' 

(Bl) 

12088 3387 772 272 624 287 184 1192 647 752 432 
"lil 2]~^^ VSb^^T^^ 75^ " 245ai ^ 35ai ^ 1575ai ^ 1225a7' 

(B2) 



where the 7^ parameters are defined in Eq. (^). For a = 1, these expressions reduce to the initial kurtosis 
factors given by the ZA, and at late times they approach the asymptotic exact values (Bernardeau 1994). 
Similarly, for EDFW initial velocities we obtain instead 



272 

S'4(a) = 54(00) - -— 
3o a 

r4(a = r4(oo - — 
6b a 



28 7 8 
T57 ^ 225 a2 
28 7 8 
TSa ^ 225 a2 



5248 32 7 



735 02 
4768 

7" 

735 a2 



15 a2 
256 7 



5056 

4725 at 
2528 



1024 
3675 a^' 
768 



105 ai 1575 ai 1225 a^' 



(B3) 
(B4) 



where 64(00) and T4(oo) denote the asymptotic exact-dynamics values given by the first four terms in 
Eqs. (Bl). For p = 5, Eq. (pS) yields for the pentosis parameters from standard ZA initial velocities 



200575880 1847200 7 6940 7^ 235 7^ 1490 72 50 7 72 10 73 
305613 3969 ^ 63 27 63~ ^ 9 27^ 
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94262120 161440 7 260 7^ 235 7^^ 130 72 50 7 72 10 73 
^^^^ ~ 305613 567 3 27 7 9 W 

14792 4108 7 94 7^ 20 72 5512 1352 7 1384 68888 1012 7 26 7^ 

49 a 21a So"''' 3 a "''63a2~ 45 a2 ~ 225 a^ ~ 1029 ai ^ 21 ai 3ai 

40 72 92072 26248 7 8752 1944 96 7 1088 2592 

21 3969 ai 2835 ai 5775 a-r 343 35 a^ 1225 a^ 8575 a^' 

whereas for EDFW initial velocities we get 
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where ^5(00) and T5(oo) denote denote the asymptotic exact-dynamics values given by the first hne in 



Eqs. (B5) and (BC) 



C. Results for Transients from Second- Order Lagrangian PT Initial Conditions 

To calculate the properties of initial conditions in 2LPT, it is convenient to take advantage of the results 
by Munshi, Sahni & Starobinsky (1994), who derived the density and velocity divergence vertex generating 
functions in 2LPT for the case in which smoothing effects are neglected. The effects of smoothing can then 
be included by the mapping given by Bernardeau (1994b). We refer the reader to these papers for details, 
here we just present the results of these calculations including the effects of transients. For simplicity, we just 
display results assuming 7^ = for p > 2. The skewness parameters, ^3 and T3 do not show any transient 
behavior in 2LPT, since these are reproduced exactly. For the kurtosis factors we have: 

54(a) 
Tiia) 

whereas for p — 5 we obtain: 
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(C5) 



(C6) 



Finally, the case p — 6 yields: 

16275904 134872 7 7544 7^ 108800 4160 7 50888 5567104 
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D. Lagrangian Perturbation Theory 

In this appendix we present the results needed to implement 2LPT initial conditions in numerical 
simulations. These perturbative results are not new and have been reported before in the literature (e.g. 
Buchert et al. 1994, Bouchet et al. 1995), but they are collected and summarized here with emphasis on 
the practical issues to ease the implementation by interested readers. In particular, Buchert et al. (1994) 
present the necessary results to implement third-order Lagrangian PT (3LPT) solutions. However, a similar 
calculation to that in Appendix ^ shows that this only improves over 2LPT by a factor of two or so in the 
expansion required to erase transients. Given the additional complexity of 3LPT, which involves solving 
three additional Poisson equations, it does not appear worthwhile to consider. 



D.l. Basic Results of Second-Order Lagrangian PT 

In Lagrangian PT, the object of interest is the displacement field ^{q) which maps the initial particle 
positions q into the final Eulerian particle positions x, 

x = q+^{q). (Dl) 
The equation of motion for particle trajectories x{t) is 



(Px ^ ^ dx 

where $ denotes the gravitational potential, and V the gradient operator in Eulerian coordinates x. Taking 
the divergence of this equation we obtain a closed equation for the displacement field 

where we have used Poisson equation together with the fact that 1 + 5{x) — J^^, and the Jacobian J{q,T) 
is the determinant 

J(q,r) =Det (5,, (D4) 
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where '^i j = d^i/dqy Equation ( p^ ) can be fully rewritten in terms of Lagrangian coordinates by using 
that Vi — {Sij + ^i,j)~^'Vq^, where Vq = d/dq denotes the gradient operator in Lagrangian coordinates. 
The resulting non-linear equation for ^(q) is then solved perturbatively, expanding about its linear solution, 
the Zel'dovich (1970) approximation 



Vg-*W =-A(r) %), (D5) 

which incorporates the kinematic aspect of the collapse of fluid elements in Lagrangian space. Here 6{q) 
denotes the (Gaussian) density field imposed by the initial conditions and Di{t) is the linear growth factor, 
which obeys Eq. (0). The solution to second order describes the correction to the ZA displacement due to 
gravitational tidal effects and reads 



V, • = ii?2(r)^(*[>g - *SfS), (D6) 

(e.g., Bouchet et al. 1995) where D2{t) denotes the second-order growth factor, which for 0.1 < ft < 3 
(A = 0) obeys 



D,{r)»^^Dl{r)n-y'^^^^Dl{r) (D7) 

to better than 0.5% and 7%, respectively (Bouchet et al. 1995), whereas for flat models with non-zero 
cosmological constant A we have for 0.01 < fl < 1 

D2{r) « -^^Dlir) f^-^ua ^ ^^_Dl{r), (D8) 

to better than 0.6% and 2.6%, respectively (Bouchet et al. 1995). Since Lagrangian solutions up to second- 
order are curl-free, it is convenient to define Lagrangian potentials 0^^-' and 0^^^ so that in 2LPT 

x{q)^q-Di V,0(i) + D2 V,^^^) ^ ^g) 
and the velocity field then reads {t denotes cosmic time) 

where H is the Hubble constant, and the logarithmic derivatives of the growth factors /,; = (d In 1?^)/ (din a) 
can be approximated for open models with 0.1 < < 1 by 



/i ~ f^'/', h ~ 2 (Dii) 

to better than 2% (Peebles 1976) and 5% (Bouchet et al. 1995), respectively. For flat models with non-zero 
cosmological constant A we have for 0.01 < < 1 



/2«2f7^/", (D12) 

to better than 10% and 12%, respectively (Bouchet et al. 1995). The accuracy of these two fits improves 
significantly for > 0.1, in the range relevant for the present purposes. The time-independent potentials in 
Eqs. (D9) and (DIO) obey the following Poisson equations (Buchert et al. 1994) 
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%4>''Kq) = S{q), (DlSa) 

v^</,(2)(q) = E['^S^(9)0S(g)-('^g(g))'], (D13b) 



D.2. Setting up Second-Order Lagrangian PT Initial Conditions 

We now have all the elements to describe a simple prescription to implement 2LPT initial conditions 
(hereafter, a tilde denotes a Fourier-space quantity): 

- A Gaussian density field, S(k) is generated in Fourier space from the desired linear power spectrum, 
and therefore ^'-'^^(fc) follows from Eq. (DlSa) in Fourier space. 



- Inverse fast Fourier transform (EFT) of 0(^)(fe) yields (/''■^-'(g) on the grid of unperturbed particle 
positions. 

- Vg0^^^ is obtained by differencing (^'■^•'(q) along the three directions, giving the necessary ingredients 
to displace particles and assign velocities according to the ZA. So far, this is the usual procedure. 

- The array Vqt/''-^-' can be stored temporarily in the array reserved for the velocities assignment. Then 
we need three additional arrays of dimension iVgrid (where A^grid is the dimension of the grid used to set up 



the initial conditions, usually equal to the number of particles) to calculate the source term in Eq. (D13b). 
By differencing the components of the array V^^^^-* in a diagonal fashion we obtain the diagonal terms 
V|2</'^^\ and V§3(/)^^\ each stored in a Ng^d array, which are then multipHed together to give the 



first contribution to the source term in Eq. (D13b). The second contribution, consisting of the non-diagonal 



terms (f>^lj{q) in Eq. ( D13h| ), is obtained by simply differencing and accumulating in turn, without the need 



of additional memory space. This yields the source term to the second Poisson equation on the grid of 
unperturbed particle positions, stored in an array of dimension A^grid- 



- EFT of this source term, solving Eq. (DlSt) for ^^^^(fc) in Fourier space, and inverse EFT yields the 
second-order potential i/)'^^^ (q) on the grid of unperturbed particle positions. 

- The 4>^'^\q) array is then differenced along one direction (recycling one of the two usable A^giid arrays 
to thus store Vi0^^-') and then particles are displaced and velocities assigned along the direction in question 



by combining Vic/)*^^^ and Vit/)'-^-' according to Eqs. (D9) and (DlC). The same procedure is applied in turn 
to the two remaining directions. 

We therefore see that the additional requirements for 2LPT initial conditions generation over the stan- 
dard ZA scheme is basically two EFT's, a 3 x Ag^id array, and differencing of the Lagrangian potentials. 
These are indeed very modest by speed and memory considerations, especially given the fact that setting up 
initial conditions is always a small fraction of the cost of running the numerical simulation. In view of the 
substantial improvement (more than one order of magnitude) on the required amount of expansion to erase 
transients from initial conditions, it seems well worth the minimal extra effort. 
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Fig. 1. — The ratio of the tree- level Sp parameters at scale factor a to their asymptotic exact dynamics value 
for scale- free initial spectra with different spectral indices: S3 (solid lines), ^4 (dotted), S5 (short-dashed), 
Sq (long-dashed), ^7 (dot-short-dashed), and Sg (dot-long-dashed). The values at a = 1 represent those set 
by the ZA initial conditions. For cosmological models with i7 < 1, the horizontal axis becomes the linear 
growth factor. 
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Fig. 2. — Same as Figure 1, but for the Tp parameters. 
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Fig. 5. — Symbols show the ratio of the Sp parameters at scale factor a measured in SCDM numerical 
simulations (Baugh, Gaztafiaga & Efstathiou 1995) to their asymptotic tree-level exact dynamics value as 
a function of smoothing scale R. Symbols represent a = 1 (open triangles), a = 1.66 (filled triangles), 
a = 2.75 (open squares) and a = 4.2 (filled squares). Error bars denote the variance of measurements in 10 
realizations. Solid lines correspond to the predictions of transients in tree-level PT, expected to be valid at 
large scales. 
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Fig. 6. — The top left panel shows the 7p parameters as a function of smoothing scale for a F = 0.5 CDM 
model. The solid line corresponds to ncs{R) = j(R) — 3, and 72 (dotted line), 73 (short-dashed), 74 (long- 
dashed), 75 (dot-dashed). The top right panel shows the ratio of the Sp parameters to those calculated by 
setting 7p = for p > 2 as a function of scale. Line-styles denote ^4 (solid), ^5 (dotted), Sq (short-dashed) 
and (long-dashed). The left bottom panel shows the analogous calculation for a F = 0.21 CDM model. 
Right bottom panel presents a similar result for the Tp parameters in the F = 0.5 CDM model. 
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Fig. 7. — The ratio of the tree-level Sp parameters at scale factor a to their asymptotic exact dynamics 
values for CDM models as a function of scale factor a. Top panels correspond to a F = 0.5 CDM model at 
smoothing scales R = 10, 100 h~^Mpc, and bottom panels show the corresponding results for a F = 0.21 
CDM model. Line styles correspond to S3 (solid), S4 (dotted), S5 (short-dashed), and Se (long-dashed). 
The upper axes in the bottom panels denote the scale factor a for Cl(a) = 0.3, with A = 0. 
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Fig. 8. — Same as Fig. ^ but for the Tp parameters. 
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Fig. 9. — Same as Fig. but for initial velocities set as in EDFW, see Eq. ([l8|). The additional dot-dashed 
Hne denotes 5*7. 
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Fig. 10. — Same as||, but for initial velocities set as in EDFW, see Eq. (|l8|). The additional dot-dashed line 
denotes T7. 
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Fig. 11. — Same as Fig. ^ but for initial conditions set using second-order Lagrangian perturbation theory. 
Compare with Figs. ^ and ^ (note the difference in plot scales). 
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Fig. 12. — Same as Fig. |[ but for the Tp parameters. Compare with Figs. ^ and | 



